function PV_dgcomp(grid,base,U,k,l,outfile,test_name,varargin)
% PV_dgcomp(grid,base,U,k,l,outfile,test_name)
% PV_dgcomp(grid,base,U,k,l,outfile,test_name,diagnostics)
% PV_dgcomp(grid,base,U,k,l,outfile,test_name,diagnostics,varnames)
%
% Write a DG solution U of the compressible Navier-Stokes problem in
% a format readable by paraview (.vtu). U is a three index array: the
% first index corresponds to the variable, the second one to the
% element degree of freedom and the third one indicates the element.
%
% All the variables are interpolated on the same grid using for each
% element a representation of degree (k,l), see BuildResampleMesh.m
% for further details. The paraview output is structured as a
% discontinuous space.
%
% The optional argument varnames can be used to specify the names of
% the fields. If it is present, it is indexed with {i} and used to get
% the names of the fields provided in U and (possibly) in diagnostics.
%
% Note: paraview only understands 3D data, and the usual way to work
% with 2D data is by setting the third coordinate to 0. This function
% handles the 2D and 3D case, adding zeros as required for the 2D
% case.
%
% See also: interpolate.

 % Interpolate to degree k
 initialize_Uk = 1;
 for i=1:size(U,1)
   [Xi,Ti,Ui,cell_type] = ...
     el_resample(grid,base,permute(U(i,:,:),[2 3 1]),k,l);
   if(initialize_Uk)
     Uk = zeros(size(U,1),size(Ui,1),size(Ui,2));
     initialize_Uk = 0;
   end
   Uk(i,:,:) = permute(Ui,[3 1 2]);
 end

 % Set the variable names
 if(nargin>8)
   varnames = varargin{2};
   for i=length(varnames)+1:size(Uk,1)-grid.d+1
     varnames{i} = ['unknown ',num2str(i)];
   end
   if(nargin>7)
     for i=length(varnames)+1:size(Uk,1)-grid.d+1+size(varargin{1},1)
       varnames{i} = ['unknown ',num2str(i)];
     end
   end
 else
   % use default names
   varnames = {'density','energy','momentum'};
   for i=2+grid.d+1:size(Uk,1)
     varnames{end+1} = ['tracer ',num2str(i-2-grid.d)];
   end
   if(nargin>7)
     for i=1:size(varargin{1},1)
       varnames{end+1} = ['diag ',num2str(i)];
     end
   end
 end
 
 % Write the output file
 fid = fopen(outfile,'w');
 
 line = '<?xml version="1.0"?>';
 fprintf(fid,'%s\n',line);
 line = '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="BigEndian">';
 fprintf(fid,'%s\n',line);
 line = ' <UnstructuredGrid>';
 fprintf(fid,'%s\n',line);

 % Write the grid
 npoints = size(Xi,2)*size(Xi,3); % total number of points
 n_v4e   = size(Xi,2);    % number of vertexes per element
 n_v4c   = size(Ti,1);    % number of vertexes per cell
 n_c4e   = size(Ti,2);    % number of cells per element
 n_cells = n_c4e*grid.ne; % total number of cells
 line = [ '  <Piece NumberOfPoints="' , num2str(npoints,'%1d') , ...
                 '" NumberOfCells="'  , num2str(n_cells,'%1d') , '">'];
 fprintf(fid,'%s\n',line);

 line = '   <Points>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Float32" NumberOfComponents="3" Format="ascii">';
 fprintf(fid,'%s\n',line);
 if(grid.d==2) % add z to the data
   Xi(3,:,:) = 0;
 endif
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e ',Xi(:,i,j));
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Points>';
 fprintf(fid,'%s\n',line);

 line = '   <Cells>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="connectivity" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ne
   i_shift = (i-1)*n_v4e;
   for j=1:n_c4e
     fprintf(fid,' %i',i_shift+Ti(:,j)-1);
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="offsets" Format="ascii">';
 fprintf(fid,'%s\n',line);
 i_shift = 0;
 for i=1:grid.ne
   for j=1:n_c4e
     i_shift = i_shift + n_v4c;
     fprintf(fid,' %i',i_shift);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="types" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ne
   for j=1:n_c4e
     fprintf(fid,' %i',cell_type);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Cells>';
 fprintf(fid,'%s\n',line);

 % Write the data
 line = '   <PointData Scalars="scalars">';
 fprintf(fid,'%s\n',line);

 % density
 in = 1;
 line = ['    <DataArray type="Float32" Name="',varnames{in}, ...
         '" Format="ascii">'];
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e\n',Uk(1,i,j));
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);

 % energy
 in = in+1;
 line = ['    <DataArray type="Float32" Name="',varnames{in}, ...
         '" Format="ascii">'];
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e\n',Uk(2,i,j));
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
  
 % momentum
 in = in+1;
 line = ['    <DataArray type="Float32" Name="',varnames{in}, ...
         '" NumberOfComponents="3" Format="ascii">'];
 fprintf(fid,'%s\n',line);
 if(grid.d==2) % add Uz to the data
   for i=size(Uk,1):-1:5 % move tracer data
     Uk(i+1,:,:) = Uk(i,:,:);
   end
   Uk(5,:,:) = 0;
 endif
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e ',Uk(3:5,i,j));
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);

 % Tracers
 for ii=6:size(Uk,1)
   in = in+1;
   line = ['    <DataArray type="Float32" Name="',varnames{in}, ...
           '" Format="ascii">'];
   fprintf(fid,'%s\n',line);
   for j=1:size(Xi,3)
     for i=1:size(Xi,2)
       fprintf(fid,'%.15e\n',Uk(ii,i,j));
     end
   end
   line = '    </DataArray>';
   fprintf(fid,'%s\n',line);
 end

 % Additional diagnostics
 if(nargin>7)
   D = varargin{1};
   for i=1:size(D,1)

     [Xi,Ti,Ui,cell_type] = ...
       el_resample(grid,base,permute(D(i,:,:),[2 3 1]),k,l);
     Uk(1,:,:) = permute(Ui,[3 1 2]);

     in = in+1;
     line = ['    <DataArray type="Float32" Name="',varnames{in}, ...
             '" Format="ascii">'];
     fprintf(fid,'%s\n',line);
     for j=1:size(Xi,3)
       for ii=1:size(Xi,2)
         fprintf(fid,'%.15e\n',Uk(1,ii,j));
       end
     end
     line = '    </DataArray>';
     fprintf(fid,'%s\n',line);

   end
 end


 line = '   </PointData>';
 fprintf(fid,'%s\n',line);

 % Close the grid
 line = '  </Piece>';
 fprintf(fid,'%s\n',line);

 % Close the file
 line = ' </UnstructuredGrid>';
 fprintf(fid,'%s\n',line);
 line = '</VTKFile>';
 fprintf(fid,'%s\n',line);

 fclose(fid);
return

